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Abstract. An early-stage version of simulation package 'ELSES' (Extra-Large- 
Scale Electronic-Structure calculation) is developed for electronic structure 
and dynamics of large systems, particularly, nm-scale or lOnm-scale systems 
(www.elses.jp). Input and output files are written in the Extensible Markup 
Language (XML) style for general users. Related pre-/post-simulation tools 
are also available. Practical work flow and example are described. A test 
calculation of GaAs bulk system is shown to demonstrate that the present code 
can handle systems with more than one atom species. Several future aspects are 
also discussed. 
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1. Introduction 

Process (molecular-dynamics, MD) calculation with electronic structure is essential 
as analysis or prediction tools of nanomaterials, particularly, materials in nanometer 
or ten-nanometer scales. Since structure and function of materials in these scales 
are determined by the competition among different regions, typically surface and 
bulk (inner) regions, the theory should reproduce such a nanoscale competition by 
describing correct electronic structures at different regions. For years, we have 
developed a set of theories and program code for such nanoscience researches. 
[U [HI [31 HI El El El [HI One crucial point is that large-scale quantum-mechanical 
calculation can be realized, in principle, by calculating the one-body density matrix, 
instead of one-electron eigen states, since the computational cost can be drastically 
reduced. [10] An overview of these theories can be found in the introduction part of 
Ref. [7] . Practical methods were constructed as solver methods of the one-body density 
matrix or the Green's function for a given Hamiltonian matrix. We note that some 
of the theories are purely mathematical ones, iterative linear-algebraic algorithms for 
large matrices and, therefore, should be useful in other fields of physics. Actually, 
one method, called 'shifted conjugate-orthogonal conjugate-gradient method', [6] was 
applied to an extended Hubbard model for L^-zSr^r^C^. [IT] Another crucial 
point is to construct algorithms for efficient parallel computations. Since multi- 
core CPU architectures are now built in standard workstation or personal computer, 
parallel computations are essential for actually all the computational systems. The 
calculations are realized with Slater-Koster-form (tight-binding) Hamiltonians and test 
calculation was carried out with 10 2 -1 7 atoms with or without parallelism. See FigQ] 
(a) for a bench mark, in which the parallelism was realized by the OpenMP directive 
(www.openmp.org). As a benchmark with a recent multi-core CPU architecture, we 
have tested our code with a standard workstation with four dual-core CPU's (Opteron 
2GHz), for liquid carbon with 1728 atoms. We adopted a typical Hamiltonian of 
carbon system. |f 2] As a result, a computational time is six seconds per time step in 
the process (MD) calculations and a parallel efficiency is more than 90%. Electronic 
property, such as density of states, is also calculated. [7 

Now the code has named Extra Large Scale Electronic Structure calculation 
(ELSES) code (www.elses.jp). It is reorganized as a simulation package 
with input/output files in the Extensible Markup Language (XML) style 
( |http:/ /www. w3.org/XML/), such as shown in FigQ] (b), for a wider range of users 
and applications. This article describes structure, example and future aspect of the 
simulation code . 

2. Structure of code 

Structure of the code is illustrated schematically in Fig. [2l The outer loop of the code 
is the loop of time evolution in MD simulation (MD part) according to the Newton 
equation; 

MjRj = Fj (1) 

where Mj,Ri are mass and position of 7-th atom (ion), respectively and Fj is the 
force on the atom. The time evolution is performed numerically by a finite-difference 
method, as in classical MD simulations. At each time step, the quantum-mechanics 
(QM) part is called from the MD part, in which the force is calculated from the 
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Figure 1. (a) The computational time as a function of the number of atoms (N). 
El El E] The time was measured for metallic (fee Cu and liquid C) and insulating 
(bulk Si) systems with up to 11,315,021 atoms, by the conventional eigenstate 
calculation (EIG) and by our methods for large systems; Krylov-subspace method 
with subspace diagonalization procedure (KR-SD), and Wannier-state method 
with variational and perturbative procedures (WS-VR, WS-PT). See the original 
papers [2] \E\ [7] for the details of parallel computation, (b) An example of XML 
formatted input file. In this example of input file, calculation conditions are set 
as the initial atomic structure file: C60.xml, integration algorithm for molecular 
dynamic calculation: velocity Verlet, time step: l.Ofsec, total simulated time: 
250fsec. Results of atomic structures are written in a file named 'C60-position.xyz 
' every 4 steps, and restart file 'C60-restart.xmr is updated every 10 steps. 
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Figure 2. Illustrated structure of the code 



position of atoms; 

=> H => p => {Fj}. (2) 

Here H is Hamiltonian matrix and p is the one-body density matrix, defined formally 

as 

p(ry)^M*(r)<Mr'), (3) 

i 

from eigen states {4>i(r)} and their occupation number {/»}. The diagonal elements 
(r = r') gives electron density (n(r) = p(r,r)) and the off-diagonal ones (r ^ r') are 
responsible for quantum-mechanical effect. In an order- N calculation, such as Krylov- 
subspace method [4J [7], the density matrix is calculated without eigen states {<f>i(r)}. 
In the present article, the calculation method of the density matrix from Hamiltonian 
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(H =>■ p) is called 'solver'. In the practical code, the matrices are given within atomic 
orbital representation 

H{ot,f3,I,J) — (a,I\H\jp) (4) 
p(a,/3,I,J) = (a,I\p\J/3) (5) 

where the suffices /, J denote atoms and the suffices a, j3 denote orbitals. The local 
density of states (LDOS) is calculated through the Green function G(e), with a Krylov- 
subspace solver [6]; 

H^G^ LDOS, (6) 

which is carried out in a post-simulation tool. 



- Atomic structure in non-XML-style file 
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Figure 3. Schematic figure of the work flow of research with examples of results 
(right panels). As examples, silicon cleavage process, [5] formation process [8] of 
helical multishell gold nanowire | 13| . local density of state of liquid carbon [7] 
appear, at the upper, middle and lower panels of the right side, respectively. 

Figure [3] indicates the work flow in research; Two input files in the XML style 
are needed; (a) file calculation conditions (See Fig.QJb)) and (b) file for initial atomic 
structure. Several conversion tools, as pre-simulation tools, are available so as to 
create the XML-style file of atomic structure from non-XML style files, such as files 
in the conventional XYZ style. 

The XML style is used commonly among text-based electronic files for sharing 
data through the internet. It is flexible and extensive, because it allows us to define 
their own items. In an XML-style file, each item is designated by start tag (<•••>) 
and end tag (</•••>). As an example of our XML-style file, Fig. [1Jb) contains a 
line of '<temperature unit="kelvin">273.15 </temperature>', which means that the 
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temperature of the system is set to be T=273.15 Kelvin. Other units can be used by 
rewriting the 'unit' part, such as l unit="eV" ', for eV. 

We note that the implementation of a XML-style input file is important in 
practical simulations of nanostructure materials, since various conditions are required. 
For example, the fracture simulation of silicon nanocrystal [5] was realized by imposing 
an external load on the atoms in a limited region near the sample boundary. The 
extendibility of XML-style file can satisfy these detailed conditions, by adding newly- 
defined tags for its own purpose. 

After the MD simulation, we analyze the atomic processes and electronic structure 
in detail as well as to visualize atomic structures. The LDOS calculation tool 
is available as a post-simulation tool, as discussed above. The tools should be 
developed further, particularly, for analyzing electronic structure, such as crystal 
orbital Hamiltonian populations |14[ [6], a quantitative visualization method of a 
chemical bond from the energetics with off-site elements of the Green's function G(e). 

3. Test calculation of compound; example of bulk GaAs 

Bulk GaAs was calculated so as to demonstrate that the present code can handle 
systems with more than one atom specie. We adopt a Slater-Koster-form Hamiltonian 
of GaAs with s, p and s* atomic orbitals. |T5] The atomic energy level of the s* orbital 
is located within the conduction band and its physical origin is a spherical average 
of the five d orbitals. The formulation of s, p, and s* orbitals was introduced [16] 
among various semiconductors, for reproducing the valence band and the bottom of the 
conduction band and was used in papers, such as Refs. [T51 El HH1 HH HOI 123 12U > for 
liquid, amorphous, defect, surface and quantum dot. Figure 4 shows calculation results 
of bulk GaAs, in which the cubic periodic cell with 64 atoms is used. Here the Krylov- 
subspace method with subspace diagonalization [H [Jj is adopted for solver routine of 
the density matrix. The dimension of the Krylov subspace (Krylov dimension) should 
be set as a controlling parameter that determines accuracy and computational cost. 
The computational time is proportional to the Krylov dimension and the calculation 
will be converged to the exact one, when the Krylov dimension increases. See Ref. [7J 
for detail. Figure 4 plots the optimized lattice constant and the energy as the function 
of the Krylov dimension. Figure 4 indicates that the calculation is well converged with 
the Krylov dimension of 30; the deviations in the lattice constant and the energy are 
less than 0.01 % and less than lmeV per atom, respectively. We note that an excellent 
convergence at the Krylov dimension of 30 was found in the other systems. [H [7J 

4. Summary 

Large-scale electronic structure calculation code is being developed as a simulation 
package with the name of ELSES (www.elses.jp). For a better user interface of our 
simulation code, we have created the input/output interfaces of XML-style files. The 
pre- and post-processing tools have been also prepared for modeling and detailed 
analysis of the atomic structures. We have also confirmed that the present code can 
handle system with more than one atom specie by calculating bulk GaAs. Although 
the present stage of the simulation package is still in an early one, we believe that 
the code will provide fruitful simulations for researchers in nano-material science. As 
a future aspect, the method is being extended by implementation of general Slater- 
Koster-form Hamiltonians for wider range of materials, in which an explicit charge 
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Figure 4. The optimized lattice constant and the energy for bulk GaAs as the 
function of the Krylov dimension. 



selfconsistent treatment |23j is included. Non-equilibrium current and other electronic 
properties arc also crucial for nanoscience and should be investigated with the present 
methodologies. 
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